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A new scheme for computing dynamical evolutions and gravitational radiations for intermediate- 
mass-ratio inspirals (IMRIs) based on an effective one-body (EOB) dynamics plus Teukolsky per- 
turbation theory is built in this paper. In the EOB framework, the dynamics essentially affects the 
resulted gravitational waveform for binary compact star system. This dynamics includes two parts. 
One is the conservative part which comes from effective one-body reduction. The other part is the 
gravitational back reaction which contributes to the shrinking process of the inspiral of binary com- 
pact star system. Previous works used analytical waveform to construct this back reaction term. 
Since the analytical form is based on post-Newtonian expansion, the consistency of this term is 
(3JT)' always checked by numerical energy flux. Here we directly use numerical energy flux by solving the 

[""t \ Teukolsky equation via the frequency-domain method to construct this back reaction term. And the 

, conservative correction to the leading order terms in mass-ratio is included in the deformed-Kerr 

metric and the EOB Hamiltonian. We try to use this method to simulate not only quasi-circular 
adiabatic inspiral but also the nonadiabatic plunge phase. For several different spinning black holes, 
we demonstrate and compare the resulted dynamical evolutions and gravitational waveforms. 
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I. INTRODUCTION 

> . 

Intermediate-mass-ratio inspirals (IMRIs): compact objects with stellar mass spiralling into intermediate mass 
(50 < M < 1000M Q ) black holes are expected to be a key source for LIGO, VIRGO and the future Einstein 
telescope. It was estimated that Advanced LIGO can sec 3-30 such events per year out to distances of several hundred 
Mpc[l|. Due to a small mass ratio of IMRIs (/i/M ~ 10~ 2 — 10~ 3 ), the small body lingers in the large black holes 
strong-curvature region for about hundreds wave cycles before merger, and then produce a good signal-to-noise ratio. 
Considering the waveform-match technology of LIGO's data analysis, ones should first model enough accurate and 
numerous waveform templates for IMRI sources. 

The post-Newtonian approximation has completely broken down in this high relativistic regime. Up to now, 
numerical relativity still encounters some difficulties to simulate the binary black hole systems when the mass ratio 
becomes small. The newest result achieves the mass ratio 0.01 (without spin) only with two orbit's evolution Q. 
For extreme-mass-ratio inspirals (EMRIs), the mass ratio about 10~ 4 — 10~ 7 , the small objects can be treated as a 
5^ \ small perturbation to the large black holes gravitational field and can be modeled by black-hole perturbation theory 
built by Regge, Wheeler and Zerilli in Schwarzschild spacetimed, H| and by Teukolsky in Kerr background [f| |(|. For 
example, a successful approach is that the orbit of small body is just geodesies of test particle in the background of 
central black hole, and the spiral of small object is simulated as a sequence of adiabatically shrinking geodesies by 
changing the constants of motion due to gravitational radiation 043 • But, when mass ratio increases, for IMRIs, the 
small object cannot be treated as test particle and the self-gravitational potential of small body becomes considerable. 
And perhaps the linear perturbation approximation is not valid in this case. 

In last decade, using effective one-body (EOB) resummed analytical formalism [r7H23| to develop accurate analytical 
models of dynamics and waveforms of black hole binaries has gotten some exciting results |24l433 |. Recently, some 
researches used EOB to model the dynamical evolution of IMRIs and extract waveforms via a multipolar Regger- 
Wheeler- Zerilli type perturbation approach [33|-|35j. And an effort of constructing the EOB waveform models for 
EMRIs by calibration of several high-order post-Newtonian parameters was also presented [HI, H3] ■ 
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The framework of effective one-body has two building blocks which are dynamics and the corresponding waveform. 
But the dynamics essentially affects the resulted gravitational waveform for binary compact star system. This dy- 
namics includes two parts. One is the conservative part which comes from effective one-body reduction. The other 
part is the gravitational back reaction which contributes to the shrinking process of the inspiral of binary compact 
star system. Previous works used analytical waveform to construct this back reaction term. Since the analytical form 
is based on some assumptions, the consistency of this term is always checked by numerical energy flux. Differently 
we directly use numerical energy flux to construct this back reaction term in this work. 

Considering that the perturbation method works well for intermediate mass ratio binaries (38i l39j . in this paper, 
we model the gravitational wave sources of IMRIs by combining the EOB frame with the Teukolsky perturbation- 
based energy fluxes. Thanks to the EOB theory, we treat the binary system by an effective particle mass fi = 
mim2/(m 1 -I- m 2 ) and a deformed Kerr (or Schwarzschild if without spin) spacetime with mass M = m\ + m 2 where 
mi is the mass of massive black holes and m 2 the mass of small body. The gravitational energy fluxes and waveforms 
are calculated by solving the Teukolsky equation in frequency domain not by the EOB fluxes and waveforms, as well 
as the inspiral dynamics is computed not following the adiabatic approximation but by the EOB dynamics which was 
developed in Refs.(33j-[35j. 

We must emphasize here that the scheme we used in this paper is different with the works by Yunes et al. which 
used the Teukolsky fluxes to calibrate EOB waveforms [36i l37j . and also different with the references (33l435| which use 
perturbation method to read the waveforms but not to radiation reaction. Our calculation covers both quasi-circular 
inspiral and nonadiabatic plunge without merger and ringdown stages. The results in this paper should be important 
for checking the EOB PN-expanded energy flux. With our method, one can obtain information for IMRIs which is 
not yet been accessible by full numerical relativity. Moreover, our method may be expanded to elliptic and even 
non-equatorial orbits. In the future works of this series, we will expand this calculation to solve Teukolsky equation 
in time domain and model the entire process of IMRIs including merger and ringdown. 

This paper is organized as follows. In the next section, we introduce EOB Hamiltonian and dynamics briefly. Then 
we describe our numerical technologies used in solving the Teukolsky equation to obtain fluxes and waveforms. In 
Sec. IV we present our calculation results for several IMRI cases. Finally, conclusion and discussion of this work are 
given in the last section. 

We use units G = c = 1 and the metric signature (— , +, +, +). Distance and time are measured by the effective 
Kerr black-hole mass M. 



The well-known EOB formalism was first introduced by Buonanno and Damour about ten years ago to model 
comparable-mass black hole binaries [l7l [l8j . and was also applied in small mass-ratio systems[33 35, 40]. Here we 
assume an IMRI system with central Kerr black hole mi and inspiralling object m 2 which is restricted on the equatorial 
plane of mi, total mass M = mi + m 2 , reduced mass /i = mim 2 /M and symmetric mass ratio v = fj,/M. Here we 
don't consider the spin of the small object. Then the deformed Kerr metric takes the form [32| 
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where 

E = r 2 + a 2 cos 2 6, (6) 

At =r 2 lA(u) + ^u 2 }, (7) 

A r =A tJ D- 1 (u), (8) 

A f = (r 2 + M 2 a 2 ) 2 - M 2 a 2 A t , (9) 

, f , fd aM 3 fd Ma 3 

uifd = 2Mar + uj J 1 v v j (10) 

r r 

A{u) = 1 - 2u + 2vu 3 + v(^— - ^tt 2 )" 4 , (11) 
(u) = 1 + log[l + 6i/w 2 + 2(26 - Zv)vu 3 ] . (12) 

Note that a is the so-called deformed-Kerr spin parameter, which is not exactly equal to the spin parameter of Kerr 
black hole itself. The values of oj{ d and ui ( d g iven by a preliminary comparison of EOB model with numerical relativity 
results are about -10 and 20 respectivelv[4ll |42|. 
The EOB Hamiltonian takes the form 



where u = Mj r and 



H EOB = My/1 + 2v(H cS / f j. - 1). (13) 

The effective Hamiltonian H c h in this paper is just the Hamiltonian of a nonspinning (NS) particle in the deformed- 
Kerr metric 



H NS = (3 l Pl + ayjfi 2 +j^PiPj, (14) 

and 

a =^=' ( 15 ) 
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With all of this, the effective Hamiltonian pT)) can be written as 




A t V A t J V At sin 2 0^ B 

The EOB dynamical evolution equations under radiation reaction can be given as 

A t A r p r 



A t E 2 p 



(19) 

*-^w£ sn * (20) 

p r = (dH m /dr)/H EO B, (21) 

P* = ( 22 ) 

where = i?NS ~P4>~a~- (of order O(^)) represents the non conservative radiation reaction force. Following the 
expression given in [34L l35l |43j , we have 

32 

J" = -— ^ 5 r 4 / DIN , (23) 
o 

where /din = F lmax /-Fkcwt is the Newton normalized energy flux up to multipolar order l max . In the EOB framework, 
F im ™ is given by "improved resummation" technique of Ref. |43j for a nonspinning binary and Ref. [23j for a spinning 
one. In this paper, we use energy flux calculated by the Teukolsky perturbation theory instead of EOB energy flux 
in Eq. (|2"2")l . As discussed in Ref.|34|, the EOB dynamical equations (|19H22j) do not make any adiabatic approximation 
which is extensively used to model EMRIs by solving Teukolsky equation [lol - |l4| . 
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III. TEUKOLSKY EQUATION AND NUMERICAL CALCULATION 
A. The Teukolsky equation 

The Teukolsky formalism considers perturbation on the Weyl curvature scale 1/14 instead of metric perturbation like 
the Regge- Wheeler- Zerilli equations. For Schwarzschild black holes these two formalisms are equivalent, but Teukolsky 
formalism facilitates us to study super-massive spinning black hole. Different to time domain method in [33 35], we 
decompose -04 in the frequency domain 

/+00 
duj^RimUr) -2SZ(9)e im ^e-^\ (24) 

where p = —l/(r — ia cos 6) . The function Ri mu] (r) satisfies the radial Teukolsky equation 

2 d f 1 dRi T , 



A — I ) - V(r)R lmu = -% mu (r), (25) 

where ^ mw (r) is the source term, and the potential is 

K 2 +4i(r-M)K , , s 

V(r) = — K - '— + 8iujr + X, (26) 

where K = (r 2 + a 2 )u — ma, A = E\ m + a 2 uj 2 — 2amw — 2. The spin- weighted angular function _2>S , ^(0) obeys the 
following equation, 



1 d ( d _ 2 SZ 
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The radial Teukolsky Eq. pB"]) has the general solution 
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where the R^ lu {r) and Ri mu {r) are two independent solutions of the homogeneous Teukolsky equation. They are 
chosen to be the purely ingoing wave at the horizon and purely outgoing wave at infinity, respectively, 

Rf m ^B^ 2 e-^\ r^r + 

Rf m . = BZir 3 e^*+r-'BZ^re'^*, r -> 00; (29) 

and 

^ w = C^' + A 2 At^e"^* , r r + 

R^=r 3 DZ^ r *, r^oo, (30) 

where k = uj — majlMr^- and r* is the "tortoise coordinate". The solution (|28j) must be purely ingoing at horizon 
and purely outgoing at infinity. That is, 

Ri m Ur -> 00) = Z$ nu r l e iur *, (31) 
i?/ mw (r^r+) = Z^A 2 e -^". (32) 

The complex amplitudes are defined as 

2iuB™J r+ A(r>) 2 
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The particle in this paper is in circular orbit on the equatorial plane, thus the particle's motion is described only 
as the harmonic of the frequency ft^ ([20)1 . We define 



Then Zfj°" are decomposed as 



oj m = mO^,. (34) 
Z^ u = Z^6(uj-u; m ). (35) 



The amplitudes Z^^ fully determine the energy fluxes of gravitational radiations, 



I ^H,oo |2 

£,oo,h = i i^ i ( 36 ) 
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and the gravitational waveform: 
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In the next subsection, we briefly introduce the traditional numerical technique solving the Teukolsky equation by 
Sasaki-Nakamura translation |44j. And in Sec. Ill C, we introduce a new numerical method solving the Teukolsky 
equation by Fujita and Tagoshi|45l |4o| . 



B. The Sasaki-Nakamura equation 

As mentioned in the above subsections, in order to get Z,'°° , we should integrate the homogenous version of Eq. (1251) . 
But there is a difficulty when one numerically integrates Eq. (p?5)) due to the long-range nature of the potential V(r) in 
([25]) . In order to solve this problem, Sasaki and Nakamura developed the Sasaki-Nakamura function X(r), governed 
by a short-ranged potential, to replace the Teukolsky function R(r) The Sasaki-Nakamura equation reads as 

fpip _ F{r) ^L _ U(r)X lmu = 0. (38) 
ar* z ar* 

The functions F(r), U(r) can be found in Ref.[44]. The Sasaki-Nakamura equation also admits two asymptotic 
solutions, 

X l H mu =e- i *> r % r^r + , (39) 

X^=A%£ a P(r)e iu "+Afa u P(r)e- iur * i r -> oo; (40) 

and 

= Cr:LHr)e tpr * + CZ u P{r)e-^\ r -+ r+, (41) 

X^=P(r)e- iur *, r->oo, (42) 

where P{r), P(r) can be found in Refs.[47l. [48j. 

The solution of Eq. (|38|) is transformed to the solution of the Teukolsky equation by 



PA Ax "^ P d Ax "u 



A J Vr 2 + a? A dr ^ 



2 



(43) 



The relations between the coefficients of the Sasaki-Nakamura function and the Teukolsky function are 



/tin 1 AiJ 1 

and the functions a, ft, T], and di muj are given explicitly in Ref.[44j and in Appendix B of (Toj . 

Our numerical scheme to integrate the Sasaki-Nakumura equation has been used in Ref.(49j for spinning test 
particle. The integrator we adopted is Runge-Kutta 7(8) and a variant of Richardson extrapolation is implemented 
to accurately compute the value of X lr ^ x> in Eqs. (|39H4"l?)) approaching to infinity and horizon. The calculation of 
spin- weighted spheroidal harmonics in Eq. (|27[) and more details can be founded in Ref . [l(]| . 
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TABLE I: The gravitational wave luminosity at infinity for r = 6M and a = 0.9 (The data labeled "Ref. [45l|" are kindly offered 
by Dr. Fujita). 
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|m| 


Ref. [45] 




This paper 




Numerical integration Relative error 
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4.61839129214686 x 10" 
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4.61839129214674 x 10" 
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4.6183912944 x 10" 


4 


4.98 x 10" 


10 


2 
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6.69474358662814 x 10" 


-7 


6.69474358662808 x 10" 
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6.6947435969 x 10" 


-7 


1.49 x 10 


-9 
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8.03430093731846 x 10" 


-5 


8.03430093731833 x 10" 
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8.0343009232 x 10" 


-5 


1.76 x 10" 


-9 
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2.92125227685241 x 10" 


-7 


2.92125227685261 x 10" 
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2.9212522761 x 10" 


-7 


2.40 x 10" 


10 
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4.30107466504698 x 10" 


-9 


4.30107466504728 x 10" 
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4.3010746567 x 10" 


-9 


1.86 x 10~ 


-9 
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1.72937852645569 x 10" 


-5 


1.72937852645567 x 10" 
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1.7293785258 x 10" 


-5 


4.05 x 10" 


10 
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4.01481595546757 x 10" 


-6 


4.01481595546768 x 10" 
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4.0148159516 x 10" 


-6 


9.71 x 10" 


10 
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9.65798420846838 x 10" 


-7 


9.65798420846792 x 10" 
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9.6579842187 x 10" 


-7 


1.96 x 10~ 


-9 



C. The semi-analytical method for Teukolsky equation 



Our main numerical technique in this work is a semi-analytical numerical method developed by Fujita and 
Tagoshi (45j|46Tl w hich is based on the MST (Mano, Suzuki and Takasugi) analytical solutions of homogeneous Teukol- 
sky equation[50l l5ll | . In the MST method, the homogeneous solutions of Teukolsky equation are expressed in terms of 
two kinds of series of special functions. The first one consists of series of hypergeometric functions and is convergent 
at the horizon 

Rtmu = e^(- X y s ' l ^/ 2 (l - X y^/ 2 Pin (x), (45) 
where p m is expanded in a series of hypergeometric functions as 

oo 

Pin{x) = ^2 a nF(n + V + 1 — IT, —71 — V — %T\ 1 — S — Z6 — %T\ x), (46) 

71= — OO 



and x — to{r + — r)/en, e — 2Mui,k = \/l — a 2 ,r = (e — ma)/n. The hypergeometric function F(a, /?; 7; x) can be 
found in mathematic handbooks. 

The second one consists of series of Coulomb wave functions which is convergent at infinity. The homogeneous 
solution of Teukolsky equation is 

Rc = z~ x - s {l - e^/z)- s -^ +T ^ 2 Mz), (47) 
where f v (z) is expressed in a series of Coulamb wave functions as 

U(z) = £ {-i) n ^ + ^ + S ~ ie \ n a n F n+v {-is -e,z), (48) 

— 00 v ' 

and z — u)(r — r_), (a)„ = T(a + n)/T(a), Fn(t], z) is a Coulomb wave function. Notice that both Eqs. (|46l48[) involve 
a parameter the so-called renormalized angular momentum. 

The key part of Fujita and Tagoshi's method is searching the renormalized angular momentum numerically. When 
harmonic frequency cj is under some critical values, real solutions of v can be found [45|. If uj exceeding these critical 
values, complex solutions of v can also be found [43. But when w becomes very large, it becomes very difficult to 
search the solution of v accurately. In this paper, we mainly employ the new method by Fujita and Tagoshi to solve 
the Teukolsky equation because of its high accuracy and high efficiency, but when co becomes large and difficult to 
find v accurately, we use the numerical integration method to evaluate the Teukolsky equation. 

For checking our numerical code's validity, in Table Q] and Table [TT1 we list our results by two kind of methods and 
the data published by Ref.[45|]. We find our results are almost same with Fujita and Tagoshi's, and the accuracy of 
our numerical integration of the Teukolsky equation is also good. 



IV. RELATIVISTIC DYNAMICS AND WAVEFORMS 



The mass ratio of IMRI systems we adopted is 0.01, and the central massive body is a Kerr black hole . The initial 
data (p r , pcf,) for solving Eqs. (|19H2"!?|) follows the so-called post-circular initial conditions described in [l9[ and take 
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TABLE II: The gravitational wave luminosity at horizon for r — 6M and a — 0.9 (The data labeled "Ref. [45(|" are kindly 
offered by Dr. Fujita). 



/ 


|m| 


Ref. [45] 




This paper 




Numerical integration 


Relative error 


2 


2 


-3.98206695526010 x 10" 


-(> 


-3.98206695525998 x 10" 


-ti 


-3.9820668046 x 10" 


-<; 


3.77 x 10" 


-8 


2 


1 


-7.10533640284950 x 10" 


-8 


-7.10533640284937 x 10 


-8 


-7.1053351261 x 10 


-8 


1.79 x 10" 


-7 


3 


3 


-1.17015834099274 x 10" 


-7 


-1.17015834099273 x 10" 


-7 


-1.1701585797 x 10" 


-7 


2.05 x 10" 


-7 
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-2.54494163085054 x 10" 


-9 


-2.54494163085075 x 10 


-9 


-2.5449398939 x 10 


-9 


6.68 x 10" 


-7 
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1 


-1.41731814300288 x 10" 


10 


-1.41731814300298 x 10" 


10 


-1.4173182794 x 10" 


10 


9.86 x 10" 


-8 
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4 


-4.26087026755638 x 10" 


-9 


-4.26087026755628 x 10 


-9 


-4.2608720408 x 10 


-9 


3.99 x 10" 


-7 


5 


5 


-1.66087471582563 x 10" 


10 


-1.66087471582562 x 10" 


10 


-1.6608748796 x 10" 


10 


9.64 x 10" 


-8 


6 


6 


-6.66108480342866 x 10" 


12 


-6.66108480342902 x 10" 


12 


-6.6610834965 x 10" 


12 


4.80 x 10" 


-7 



limit v — > 0, 

= L z , (49) 

Pr = -[—^]^/^ (50) 



We should point out here that usually we do not take the limit v — > (which was used in Rcfs433 35]) during the 
dynamical evolution; instead, we use the full EOB dynamics described in Eqs. (|19H2"2"l) . This makes it so we can contain 
the hrst order conservative self-force of the small body in dynamics. We list our numerical algorithm here for clarity: 

1. We calculate the quasi- circular data for the small body at time-step t^; 

2. Use the orbital frequency from the step 1 to solve the Teukolsky equation in frequency domain by the numerical 
methods we introduced in Sec. Ill up to Z max = 12; 

3. Then get the flux and waveforms (Teukolsky-based) by summing over the multipoles at the time-step tf, 

4. Calculate the dynamics of the small body for the time-step U+i by solving the EOB dynamics equations (|19M22I) ; 

5. Iterate from the step 1. 

The method of modeling back reaction and waveforms based on the Teukolsky equation in frequency domain are 
quite successful in quasi-periodic inspirals of extreme mass-ratio cases. It is a good standard for checking other 
methods. For example, in Ref. 37], Yunes et al. used the Teukolsky waveforms to construct more accurate EOB 
waveform models. But while the mass ratio fi/M ~ 10 -2 , the precision of Teukolsky energy fluxes and waveforms 
would decrease because it is first order perturbation theory and also the subleading terms in the mass-ratio introduce 
conservative correction is not included in the Teukolsky equation. It is similar in Refs. j33T-[3o| where they calculated 
/din by the EOB analytical radiation reaction in the v — > limit and gravitational waves by the Regge- Wheeler- Zerilli 
perturbation equation. Anyway, there are several reasons drive us choose the Teukolsky-based flux to replace the EOB 
analytical one to calculate the radiation reaction in this work. One is that our Teukolsky-based flux does not include 
any low velocity approximation comparing with the post-Newtonian expanded EOB flux. But we also notice that 
the /din based PN resummed waveform extends effectively the PN waveform beyond the slow velocity approximation 
quite well 52]. One is that the Teukolsky-based flux can be calculated to arbitrary multipole (but actually because 
of a computation limit, we compute to l max = 12) against the EOB's maximum I = 8 until now. The last one is that 
the black hole absorbed term is naturally included in the Teukolsky-based energy flux. 

There is also a problem to evolve the plunge orbit inside the Kerr or Schwarzschild ISCO(innermost stable circular 
orbit) with the adiabatic approximation frequency-domain code. The plunge orbit is unstable and do not have a 
well-defined frequency spectrum. But as same as the EOB-based evolution of plunge orbit, we can just use the EOB 
dynamical equation (|20[) to define the orbit frequency and calculate the harmonic frequency in Eq. (|34[) . We think this 
is an alternative way to calculate the energy flux in the Teukolsky frequency-domain frame. And the reliability of 
this definition of orbital frequency was investigated in 35] • The key point is that the orbital evolution in this paper 
is based on the EOB dynamics (I19H22[) but not the adiabatic approximation. 

First, we compare the same situation in Fig.l of Ref. [331] by our numerical algorithm with the 5PN resummed /din 
in [34|, [HI- For comparison, we take the v — > limit as same as the Eqs.(2-6) of 34]. Our results are showed in 
Figure [TJ in which both results are consistently very good at qualitative level. Because of the different methods for 
getting energy fluxes, our results of orbital dynamics and waveforms are slightly different, with the one by the 5PN 
resummed back reaction at quantitative level. In details, before ISCO, during the inspiral, both results are almost 
the same, but in plunge phase, a small but visible dephasing (about 0.16 radian) appears. Of course, this dephasing 
should be larger if involving more orbit evolution. 

Furthermore, to validate our result, we check the "convergence" of the I = m = 2 waveform modulus while v — > 0. 
The result is displayed in Fig. [5] and can clearly dedicate that the waveform we calculated has a good convergence. We 
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FIG. 1: Comparison of our method(solid line) with EOB (dashed line, the data are calculated by Dr. Bernuzzi). Left panel: 
Transition from quasi-circular inspiral orbit to plunge from ro = 7M and v = 0.01. The dashed-dotted line circle is the ISCO 
at risco = 6M and the dotted circle the horizon of black hole; Right panel: Orbital frequency fi^ versus coordinate time. 



FIG. 2: "Convergence" of the waveform when v — > for a = 0. The values of v are 0.01, 0.001, 0.0001. t — is time for the 
maximum of the waveform modulus and spatial position is the light ring at r = 3M for each value of v. 



think that this good converging trend, together with the good consistency in Fig. [T] gives a support to our theoretical 
idea which combines EOB dynamics with the Teukolsky perturbation method in the intermediate mass-ratio inspiral 
and plunge. 

Now, as an application, we show the orbits of the small object inspiralling into Kerr black holes with different spin: 
a = —0.9, 0.7. We use the full EOB dynamics equations (| 19ti22[) without the v — > approximation. This makes us 
can contain the conservative mass-ratio corrections. When a ^ 0, v ^ 0, it is no possible to write out an analytical 
expression for the angular momentum. So we use Brent's algorithm [5 3] to search L z numerically with a relative 
error less than 1CP 13 . For the Kerr black holes, because of the so-called frame-dragging effect, the radius of ISCO is 
smallcr(a > 0) or biggcr(a < 0) than the Schwarzschild case. For example, the extreme Kerr black hole, a = ±1, the 
ISCO locates at r = 1M and r = 9M, respectively. And for a — 0.7, — 0.9, nsco = 3.39313, 8.71735, respectively. 

For this reason, when a < (retrograde orbit), the small body plunges faster into the black hole. At the same 
time, when a > (prograde orbit), the small body spirals more cycles before plunge. These phenomena are clearly 
displayed in Figj3] 

In addition the waveforms are modified by different values of spin of Kerr black holes. As an example, we show 
the waveforms of a = —0.9, 0.7 cases in Figure 0] After passing the light ring, based on the analysis in Ref . (35| , the 
merge actually has started and gone into ringdown phase, the energy flux and waveform cannot be calculated by the 
frequency-domain method but can be connected by a quasi-normal mode or time-domain Teukolsky solution. For this 
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FIG. 3: Inspiralling and plunge of a small body into a Kerr black hole with mass ratio 0.01. The dashed line represents ISCO 
and dotted line the light ring. Left panel: retrograde orbit a = —0.9 beginning at r = 10M ; right panel: retrograde orbit 
a = 0.7 beginning at r — 3.5M. 




FIG. 4: GW signal shown through the I = m = 2 mode of + polarization corresponding to the dynamics depicted in Fig[3] 
The left panel is the retrograde orbit a — —0.9 and the right one retrograde orbit a = 0.7. 



reason, we stop the evolution just after the small object passing the light ring. 



V. CONCLUSIONS AND DISCUSSION 



Previous works on IMRI dynamics through the EOB method are based on approximated analytical waveforms. 
We propose an alternative method to set the back reaction term for the EOB dynamics via numerically calculated 
energy flux directly. We have presented our numerical results based on the EOB dynamics and the frequency-domain 
Teukolsky flux in this paper. This is the first try to use the frequecy-domain Teukolsky frame in the plunge phase 
of IMRIs by the help of the EOB dynamics. Our results are consistent to the previous results using the analytical 
energy fluxes of post-Newtonian approximation. 

And the results shown here are reasonable. By comparing with the full EOB evolution in deformed Schwarzschild 
space-time, a quite good coincidence before ISCO gives a support for applying the Teukolsky equation up to mass-ratio 
0.01. After ISCO, the plunge stage, a small difference comes out due to two different energy fluxes used. 

Since there are few results of post-Newtonian waveforms on elliptic, and even non-equatorial orbits, it's not easy 
to generalize the previous EOB dynamics works to these systems. In principle, our numerical energy-flux method is 
not restricted by this. The radial and polar frequencies can be obtained by evolving an "imaginary" geodesic orbit at 
every step. The "imaginary" geodesic orbit is the solution of Eqs.(19)-(22) but JF — in (22). This solution involves 
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ordinary differential equations with only few variables, so the computation time is negligible compared with the total 
evolution time. Of course, when we try to expand our method to the elliptic, and even non-equatorial orbits, the 
numerical cost will be more high. However, for EMRIs, there were several works which evolved generic orbits with 
solving the Teukolsky equation [kJ [H, EH 13 • For IMRIs, we can afford the numerical cost even for evolving a few 
tens of orbits. And more, we can generalize our analysis to generic spinning objects moving along generic elliptic, and 
even non-equatorial orbits. These would be our future works. 

The adiabatic approximation works well in inspiral with extreme mass-ratio limit, but is violated in plunge stage. 
At the same time, the EOB dynamics can be used in both inspiral and plunge evolution. And then the energy fluxes 
also are calculated by the EOB-PN expanded waveforms usually. In this traditional method, the EOB fluxes are based 
on the post-Newtonian approximation, and losses accuracy in the high relativistic region. What we do in this paper 
is, in the EOB dynamics frame, use the Teukolsky-based energy fluxes and waveforms to replace the EOB ones. In 
other words, we adopt the frequency-domain Teukolsky-based energy fluxes and waveforms but abandon the adiabatic 
approximation; we employ the EOB dynamics to evolve the binary systems but do not use the EOB energy fluxes for 
back reaction. 

This EOB plus Teukolsky frame we develop in this work can provide an alternative way to research the gravitational 
waves from IMRIs. The dynamical evolution and waveforms for several different spinning black holes are given in the 
above sections. But more detailed analysis about this EOB plus Teukolsky frame should be done farther in future, 
when more interesting results can be obtained. 
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